program mkhisqp ! -- makes histgroms from list-mode data file.
  !
  ! DESCRIPTION:
  !
  ! This is a speed-enhanced version of mkhisq. It assumes at most 8
  ! CPU cores are available for parallel computing. Since file I/O is
  ! the one that takes most CPU time, the impact of parallelization
  ! was turned out to be negligible for a ordinary PC.
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use kuty
  use suty
  use omp_lib
  implicit none
  interface
     subroutine initdecv(a,b,c)
       integer*8,intent(in),value::a,b,c
     end subroutine initdecv
  end interface
  interface
     subroutine decvq(pe,ne,pb,nb,pneo,pnbo)
       integer*8,intent(in),value::pe,ne,pb,nb,pneo,pnbo
     end subroutine decvq
  end interface
  integer,parameter::MAXPATHLEN=1024
  integer*8,parameter::nemax=5898240/4      ! maximum number of 4-byte to read
  integer*8,parameter::nlmdbufmax=737280    ! maximum number of listmode events
  integer*4,allocatable,target::lmdbufs(:,:,:)
  !
  integer,parameter::nd=8         ! total number of blocks
  integer,parameter::ndh=nd/2     ! number of blocks in a plane
  integer,parameter::mgap=4       ! number of crystals in a gap
  integer,parameter::nt=64        ! number of crystals in transaxial
  integer,parameter::na=20        ! number of crystals in axial
  integer,parameter::nx=na*ndh+mgap*(ndh-1) ! number of crystals in x
  integer,parameter::ny=nt        ! number of crystals in y
  integer,parameter::C2=ndh*(na+mgap)
  integer,parameter::nxy=nx*ny
  !
  integer*4,allocatable,target::bufs(:,:)
  integer*2,allocatable,target::hists(:,:,:,:)
  integer*2,pointer::histp(:,:),histd(:,:)
  character*(MAXPATHLEN) str,flist
  character*(MAXPATHLEN),target::fins(8),fout,foutp,foutd
  character*(MAXPATHLEN),pointer::s
  integer len
  integer*8,target::ns(8),ie,fss(8),nes(8),nbufs(8),nbs(8)
  integer,target::ndelays(8),npairs(8),istats(8)
  integer ifile,nf,ia,narg  
  integer xy1,xy2,pd

  narg=command_argument_count()
  if(narg<3) then
     call usage
     stop "*** too few arguments."
  end if
  
  fout=""
  foutp=""
  foutd=""
  flist=""
  nf=0
  nullify(s)
  
  do ia=narg-1,2,-2
     call get_command_argument(ia,str,len)
     if(len==2) then
        select case(str(2:2))
        case("o")
           s=>fout
        case("p")
           s=>foutp
        case("r")
           s=>foutd           
        case default
           call usage
           stop "*** invalid option."
        end select
        if(associated(s)) then
           call get_command_argument(ia+1,s)
           nullify(s)
        end if
        nf=ia-1
     else
        if(fout/="".or.foutp/="".or.foutd/="") exit
        call usage
        stop "*** please specify outfile correctly."           
     end if
  end do

  allocate(hists(0:nxy-1,0:nxy-1,0:1,8),stat=istats(1))
  if(istats(1)/=0) stop "*** cannot allocate memory"
  allocate(lmdbufs(4,nlmdbufmax,8),stat=istats(1))
  if(istats(1)/=0) stop "*** cannot allocate memory"
  hists=0
  ndelays=0
  npairs=0
  nbs=0
  nbufs=0
  call initdecv(na+mgap,na,C2)
  allocate(bufs(nemax,8),stat=istats(1))
  if(istats(1)/=0) stop "*** cannot allocate memory"

  !$omp parallel do
  file: do ifile=1,nf     
     call dodec(ifile)
  end do file
  !$omp end parallel do

  forall(xy1=0:nxy-1,xy2=0:nxy-1,pd=0:1)
     hists(xy1,xy2,pd,:)=sum(hists(xy1,xy2,pd,:))
  end forall

  histp => hists(:,:,1,1)
  histd => hists(:,:,0,1)

  write(*,*) "Number of files read:   ",nf
  write(*,*) "Number of bytes read:   ",sum(nbs)*4
  write(*,*) "Number of accept events:",sum(npairs)
  write(*,*) "Number of random events:",sum(ndelays)
  
  if(foutp/="") call write_hist(histp,foutp)
  if(foutd/="") call write_hist(histd,foutd)
  
  write(*,*) "Max Prompt: ",maxval(histp)
  write(*,*) "Max Delayed:",maxval(histd)
  histp=histp-histd
  where(histp<0) histp=0
  write(*,*) "Max True:   ",maxval(histp)

  if(fout/="") call write_hist(histp,fout)
 
contains

  subroutine dodec(if)
    integer,intent(in)::if
    integer tid
    integer*4,pointer::at1,at2,pd
    integer nread,nbuf_read,nlmdbuf_wrote
    integer pbuf
    integer*4,pointer::lmdbuf(:,:)
    character*(MAXPATHLEN),pointer::fin
    integer*8,pointer::ne,fs,n,nb
    integer*4,pointer::buf(:)
    integer,pointer::npair,ndelay,istat
    integer*2,pointer::hist(:,:,:)

    tid=omp_get_thread_num()+1
    lmdbuf => lmdbufs(:,:,tid)
    fin    => fins(tid)
    ne     => nes(tid)
    n      => ns(tid)
    fs     => fss(tid)
    nb     => nbs(tid)
    buf    => bufs(:,tid)
    npair  => npairs(tid)
    ndelay => ndelays(tid)
    istat  => istats(tid)
    hist => hists(:,:,:,tid)

    call get_command_argument(if,fin,len)
    fin=trim(adjustl(fin))     

    open(unit=tid+10,file=fin,status="old",access="stream",iostat=istat)
    if(istat/=0) stop "*** cannot open file."
    fs=fsize(fin)
    write(*,"(a)") "Opening file: "//trim(fin)//trim(i2a(fs))//" byte"
    if(fs==0) return
    
    ne=fs/4
    n=min(ne,nemax)     

    nread=0
    read: do        
       read(tid+10,iostat=istat) buf(1:n)
       if(istat/=0) stop "*** error during read."
       nread=nread+n
       nb=nb+n
       pbuf=loc(buf)
       dec: do
          
          call decvq(pbuf,n,loc(lmdbuf),nlmdbufmax,&
               loc(nbuf_read),loc(nlmdbuf_wrote))
          npair=npair+nlmdbuf_wrote
          sort: do ie=1,nlmdbuf_wrote
             at1=>lmdbuf(1,ie); at2=>lmdbuf(2,ie)
             pd=>lmdbuf(3,ie)
             if(pd==0) ndelay=ndelay+1
             hists(at1,at2,pd,tid)=&
                  hists(at1,at2,pd,tid)+1
          end do sort
          
          if(nbuf_read==n) exit
          pbuf=pbuf+4*nbuf_read
          n=n-nbuf_read
          
       end do dec
       
       if(nread==ne) exit          
       n=min(ne-nread,nemax)
       
    end do read
    close(tid+10)
  end subroutine dodec
  
  subroutine usage
    write(*,100)
100 format(x, "Usage: mkhis {infiles | -f flist}",/,&
         8x,"{-o outfile_true | -p outfile_prompt | -d outfile_delayed}")
  end subroutine usage
  
  subroutine write_hist(h,outfile)
    integer*2,intent(in)::h(:,:)
    character*(*),intent(in)::outfile
    integer istat
    character*128 s
    write(*,"(a,$)") "Writing "//trim(outfile)
    open(unit=1,file=outfile,access="stream",status="replace",iostat=istat)
    if(istat/=0) then
       write(*,*) "*** error during open"
    else
       write(1,iostat=istat) h
       close(1)
       if(istat/=0) then
          write(*,*) "*** error during write"
       else
          s=trim(i2a(nx))//"x"//trim(i2a(ny))
          write(*,*) "size: "//trim(s)//"x"//trim(s)//"x"//trim(i2a(sizeof(h(1,1))))//" = "&
               //trim(i2a(sizeof(h)))//" bytes"
       end if
    end if
  end subroutine write_hist

end program mkhisqp
